Probability Distribution of Residence-times of Grains in Sandpile Models 
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We show that the probabihty distribution of the residence-times of sand grains in sandpile models, 
in the scaling limit, can be expressed in terms of the survival probability of a single diffusing 
particle in a medium with absorbing boundaries and space-dependent jump rates. The scaling 
function for the probability distribution of residence times is non-universal, and depends on the 
probability distribution according to which grains are added at different sites. We determine this 
function exactly for the 1-dimensional sandpile when grains are added randomly only at the ends. 
For sandpiles with grains are added everywhere with equal probability, in any dimension and of 
arbitrary shape, we prove that, in the scaling limit, the probability that the residence time greater 
than t is exp{—t/M), where M is the average mass of the pile in the steady state. We also study 
finite-size corrections to this function. 



In 1987, Bak, Tang and Wiesenfeld (BTW) proposed 
sandpile model, as a paradigm of self-organized critical 
systems in nature 0. Since then, many different kinds 
of sandpile models have been studied. These include 
models with discrete or continuous variables 0, differ- 
ent toppling rules 0, deterministic or stochastic driving 
0, with or without local particle conservation driv- 
ing mechanisms 0, Q etc. By now, a fair amount of 
understanding of the critical steady state and critical ex- 
ponents of avalanches has been achieved, by the study of 
exactly solved models, and by numerical simulations. For 
a review of known results, see |5l^llollll|. 

However, time-dependent properties of self-organized 
critical systems have not been studied as much theoret- 
ically so far, in spite of the fact that an explanation of 
1/f noise was one of the main motivations for the ini- 
tial proposal of self-organized criticality. While a power 
spectrum of mass-fluctuations of 1 /ftype has been found 
in some 1-dimensional models |l2L |l3j , it appears that in 
higher dimensional sandpile models, the behavior is 1//^ 
hj. Piles of Ion g-grain rice have provided a very good 
experimental realization the basic ideas of self-organized 
criticality 0, The probability distribution of res- 
idence times (DRT) of grains was studied experimen- 
tally in the Oslo ricepile experiment, and by simulations 

In this paper, we will study DRT in sandpile models. 
We argue that the problem of determining DRT can be 
reduced to that of finding the probability distribution of 
hitting time of a single diffusing particle to the boundary, 
diffusing in a medium with site-dependent jump rates. 
In the scaling limit of large system sizes, DRT becomes 
a function of a single scaling variable t/L^, where t is 
the residence time, L is the linear size of the system, and 
b is some exponent. This function is non-universal, and 
is a complicated function of the spatial distribution of 
added grains used to drive the pile to its steady state. 
We determine this function explicitly for 1-dimensional 
sandpile when grains are added randomly only at the 
ends. When grains are added with equal probability ev- 
erywhere, we prove that the exact scaling function of the 



DRT is a simple exponential. This result is independent 
of dimension, and of the shape of the pile. 

Let us consider the problem in the simplest setting 
first: the BTW model m on a line of L sites, labeled by 
integers 1 to L. At each site i we have a non-negative 
integer variable Zi called the height of the pile at that 
site. The site is stable if Zi < 1. If Zi > 2, the site 
is said to be unstable, and relaxes by toppling. In this 
process, Zi decreases by 2, and Zi-i and Zi+i increase by 
1. Toppling at a boundary site causes the loss of one 
sandgrain from the pile. The pile is driven by adding 
grains at the right end, at a constant rate of one grain 
every P time steps. We assume that P is larger than the 
duration of the longest avalanche in the system, so that 
all avalanches have died before a new grain is added to 
the system. 

This model has an abelian property, and its properties 
are well- understood The long-time behavior under 
the deterministic evolution is that after an initial tran- 
sient period, it falls into a cycle of period L + 1. The 
stable configurations of the pile that belong to the cy- 
cle are L configurations having all, except one site, with 
height 1, and one configuration with all Zi = 1. If we 
start with the state with all Zi — 1, adding a particle at 
i = L gives a stable configuration with zi = 0. Adding 
a particle again, we get the recurrent configuration in 
which Z2 = 0. For each new added grain, the position 
of the zero shifts one step to the right, till after L steps 
it is at I = L. Then adding another grain, the zero dis- 
appears. We choose to say that in this case the zero is 
at i — 0. The number of toppling to get the next sta- 
ble configuration is also periodic with the same period : 
L^L-1^(L-2)...1^0^L. 

If we want to study DRT in this model, we have to 
mark the grains. However, with marked gains, the model 
is no longer abelian. This is because toppling at two ad- 
jacent unstable sites in different order no longer give the 
same result. For a full specification of the rules governing 
the motion of grains in the model, we have to define pre- 
cisely in which order the unstable sites are toppled, and 
how the grains are transferred under toppling. We choose 
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the parallel update scheme: make a list of all sites which 
are unstable at a time t, choose at random two grains 
from each of these sites ( if there are only two grains, 
both are selected), and randomly assign one of them to 
go the left neighbour, and the other to the right. All 
these grains which are to be moved are then added to 
their destined sites, at the same time. This constitutes 
a single microstep of evolution. Then we construct the 
new list of unstable sites for the next micro time-step, 
and repeat. 

The constant time elapsed between two successive ad- 
ditions of grains (P micro-steps) will be called a meso- 
time step. We measure the residence time in units of 
meso-steps. We mark all grains by the meso-time step 
when they were added to the pile. Then, if the grain 
numbered Tin (added at meso-step Tin) gets out of the 
system at meso-step Tout, we will say that its residence 
time is Tout - Tin- 

It is easily seen that the first moment of the DRT is 
the average mass of the pile. Define a variable r]{i,j) as 
the indicator function that the i-th grain is in the pile at 
the end of time meso-step j. Then, summation of ri{i,j) 
over j gives the residence time of particle i, and average 
over i gives the mean residence time. Conversely, sum 
over i, we get the mass of the pile at the end of time- step 
j , and average over j gives the mean mass of the pile |l9j . 

From our definition, it follows that the probabilities 
of different paths taken by a grain are exactly that of 
an unbiased random walker on the line. This is because 
when a grain moves under toppling, it is equally likely to 
take a step to the right, or to the left. So, for example, 
the average number of steps a grain takes before it leaves 
the pile is equal to the average number of steps a random 
walker would take from that starting point. However, 
the time between two jumps of the grain is random, and 
has very non-trivial correlations with times of previous 
jumps, and also with jump times of other particles. This 
is what makes this problem nontrivial. 

To calculate the DRT for the linear chain of L sites, 
we consider adding a marked grain at meso-time Tq. 
All other grains are unmarked, and indistinguishable. 
Then, stable configurations of the pile are in num- 
ber. Configuration in which the site a has height 0, and 
the marked grain is at site b will be denoted by Ca.b- AH 
sites other than a and b are occupied by unmarked grains. 
For each value of a, 1 < a < L, then there are L — 1 possi- 
ble configurations corresponding to different values of b. 
For the recurrent configuration with all Zi = \, we define 
a = 0, and in this case there are L possible positions of b. 
Thus there are in total possible stable configurations 
of the pile. 

Consider a particular configuration Ca.b- Adding an- 
other (unmarked) grain at i = L. li b < a, then it is 
easily seen that the wave of toppling does not reach 
the marked grain, and the final configuration is Ca+ifi- 

When b > a, the wave of toppling, started at the right 



end, reaches the site b and the site will topple. The 
marked grain will move one step to the left or right, with 
equal probabilities. If the marked grain moves to the left, 
it will move again due to toppling, unless that site has no 
grains. In this way, the marked grain can take zero, one 
or more consecutive steps to the left in one meso-step. It 
stops diffusing as soon as it takes a right step or if the 
marked grain falls on a. We thus see that on adding more 
grain, if & > a, the final configuration is Ca+i,b+Ab, with 
Ab taking values 1,0,— l...a+l — 6,a — 6 with proba- 
bilities 2-1,2-^2-3 . . . ^ 2'^-''+i, 2"-''. For (b - a) large, 
the mean square displacement ((A6)^) tends to 2. 

It is straight forward to construct the x matrix 
W giving the transition probabilities between different 
configurations. Also, knowing the probabilities of differ- 
ent configurations in the steady state, we can write down 
the probabilities of different stable configurations just af- 
ter the marked grain has been added. Then W*P(i = 0) 
gives the probabilities of different configurations at time 
t, and summing over all configurations, we get the sur- 
vival probabilities S{t) that the marked grain remains 
inside the system up to time t. In fact, using the fact 
that the position a of site with height zero changes by 
I deterministically in time, one can rewrite this problem 
in terms of (L+I) matrices Tj, j = 0, 1, 2, . . . L, where 
Tj is a (L — I) X (L — I) matrix giving transition prob- 
abilities from the L — 1 configurations with a = j to the 
L—1 configurations with a = j + 1. For a = L + 1, there 
are L stable configurations with a marked grain, but the 
configuration with b = L is transient, and hence one can 
work with the (L—1) remaining configurations. 

In an earlier paper, one of us used these to determine 
DRT exactly numerically, for L up to 150 0|. In the 
limit of large L, it was argued that Prob{t\L) tends to 
the scaling form 

Probm ^ ^fit/L'). (1) 

Here the function f{x) varies as x~^^^ for x ^ 1, and 
as exp{—Kx) for x^ 1. We now obtain the exact func- 
tional form of the scaling function /. 

Consider a marked grain at 6, at some time T, with b = 
aL,Q < a < 1, and L large. We consider the change in 
its position A6 after one cycle ( (L-|- 1) mesosteps). Fig.l 
shows the motion of grains in a cycle in one realization. 
The grain diffuses for a while, and is stuck when the zero 
is to the right of the marked grain. The average fraction 
of time it moves is a. A grain at b is hit by b waves of 
toppling 2QJ in this interval. The net displacement A6 
is sum of displacements due to these waves. Each wave 
causes a displacement with mean zero, and variance 2. 
Then by central limit theorem, the net displacement will 
be distributed normally with variance given by 26. Thus 
A6 is of order •\/2aL, and is much smaller than L when 
L is very large. Then, for times t ^ we can average 
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over the motion in a cycle, and say that if the marked 
grain is at i, it moves to the left or right neighbor with a 
rate (i/L) per unit time. If P{i, t) is the probability that 
the marked grain is at i at time the evolution equation 
for t) for times t 1 is 



(2) 

At time t = 0, we can assume that the marked particle 
is at i = L, so that P(i,t = 0) = Si^^. Integrating this 
equation, we determine the survival probability S{t) = 
P(«, t), and then the DRT is given by 



Prob{t\L) = S{t)- S{t + 1) 



(3) 



We introduce the reduced coordinate ^ — i/L, and r = 
t/L'^ and consider the Ea. HlU|) when L is large. In terms 
of these reduced variables, the evolution equation for the 
probability density P(^,t) becomes, in the continuum 
limit, 



exp{-Ct/L'^), where we get C = kl = 3.6705. This value 
is in good agreement with the value obtained by extrapo- 
lation of estimates obtained by measuring the coefficients 
of the exponential determined by exact diagonalization of 
the master equation for finite L |l9l| . 

The generalization of these results to d dimensions is 
straight-forward. We consider a d-dimcnsional sandpile 
model on a lattice with number of sites V . We assume 
that when a new grain is added, the site x is chosen 
with probability r{x). Clearly, the sum of r{x) over all 
sites is 1. In the steady state, n(a;), the average number 
of topplings at x per added grain, satisfies the equation 
(using conservation of sand grains) 



-r{x), 



(7) 



with n{x) = at the boundary. The solution of this 
equation is 



n{x) — G{x, x')r{x'), 



(8) 



(4) 



We can integrate this equation numerically using the 
initial condition i = 0) = 5{^-\ + l/L). The scaling 
function f{x) is given by. 



fix) 



(5) 



Let Lpj (^) be solution to the eigenvalue equation corre- 
sponding to eigenvalue \j 



(6) 



where (pj{^ = 1) =0 corresponding to an absorber being 
present at « = L -I- 1. At ^ = 0, we do not need to 
assume any special condition, as the absorber at i = is 
automatically taken care of by the fact that rate of jump 
out of I is 2i/L, which becomes zero at i = 0. 

We look for a solution cpj{^) that does not diverge at 
^ = 0. Expanding in a power series, and matching 

coefficients, we get 



where G{x,x') is the average number of topplings at x 
due to addition of a grain at x' , and is equal to the inverse 
of the toppling matrix A [l8| . 

The important point to realize is that while avalanches 
in sandpile can spread quite far, the typical distance trav- 
eled by one marked grain in an avalanche is much smaller 
than L. In fact, in many cases, we expect it to be of order 
1. During its motion to the boundary, the marked grain 
would be involved in a large number of avalanches. At 
time-scales much larger than a meso-step, the motion is 
diffusive, with the jump-rate out of different sites being 
space-dependent because on the average some parts of 
the lattice have more avalanche activity than others. 

Consider a grain at site x at time t. Let its position be 
x+Ax after At new grains have been added, where L'^ ^ 
At ^ 1. As the path of the grain is an unbiased random 
walk, we have ((Ax)^) = s, where s is the average number 
of jumps the grain makes in this interval. Assuming that 
\Ax\ <^ L, and that n(x) is a slowly varying function of 
X, we see that s has to be proportional to n{x)At, total 
no of toppling waves during time interval At. Let us 
say s = Kn{x)At, where K is some constant. Writing 
{{Axf) = T{x)At, we get 



where Ii {x) is the modified Bessel function of order 1 121| . 
The eigenvalues \j is obtained by imposing the condi- 
tion ip{$^ = 1) = 0. Thus if the j-th zero Ii{z) occurs 
at zL2ikj, then \j = k'j. At large times t, S{t) varies as 



T{x) = Kn{x), 



(9) 



where the constant K depends on the details of the 
model. For large times t, the probability-density P{x,t) 
satisfies the equation 



d_ 
di 



P{x,t)^-KV^[n{x)P{x,t)] 



(10) 
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with the initial condition is given by 



P{x,t = 0) = r(x). 



(11) 



It may be noted that Ea. (|10(l is not the diffusion 
equation with space-dependent diffusion constant D(x), 
where the right hand side would have been of the form 
\7[D{x)\7P{x)). The net current between two sites de- 
pends on the difference in the product nP at the two 
sites, and can be non-zero even if VP is zero. 

Solving this differential equation, with the condition 
that P{x, t) is zero at the boundary corresponding to the 
absorbing boundaries, we can determine P[x, t) at any 
time t. Integrating over x determines the probability 
that marked particle remains in the system at time 
and the DRT is obtained from the survival probability 
using Eq.jni). 

Consider, as an example, the case of a linear chain 
with L sites, when we add particles at each step at either 
of the two ends with probability 1/2. In this case, we 
get n(x) = ^, independent of x, and K — 2. One can 
then solve the Eq. lfTU|l analytically, and a straightforward 
calculation gives S{t) = 03(O,t) — 03(77/2, t), where r = 
7r^i/2L^ and 9^{z,t) is the Jacobi theta function 
defined by 

00 

^3(2, t) = 1 + 2 ^ exp(— n^T)cos(2nz). 

n=l 

In fig 2, we have plotted the analytically computed 
survival probability S{t) versus t/L"^ for this case and 
compared with the results of simulation for L = 100 us- 
ing 10^ grains. We have also shown the result of the 
numerical integration of Eq. ^ to determine the scaling 
function in the case where the grains are added only at 
one end, and compared it with the simulation data ob- 
tained for L = 100 using 10^ grains. Clearly, the agree- 
ment is excellent. This supports our arguments used to 
obtain Eq.lUHll. 

The Ea. lfTnji is very easy to solve in the special case 
when sand grains are added randomly at any sites in the 
system. Clearly here r{x) = l/V where V is the number 
of sites in the lattice. Then n{x) is a solution to the 
equation 



V'^n{x) = 1/F, 



(12) 



The function P{x,t) = T{t), for aU x, satisfies Ea. (|Tn|l in 
any dimension if 

dt 2V^' 

With the initial condition P{x, t — Q) — l/V , it's easy 
to see that the full solution is given by 



The probability of survival upto time t is VP{x,t), 
and we see that it decays in time as a simple exponential. 
Using the fact that the mean residence time is the average 
mass M in the pile, we see that K/2V — 1/M . This then 
implies that 

S{t) = exp{-t/M). (14) 

We note that our derivation depends only on Eq.Q 
and Eq.©. These two equations are valid under the 
conditions of local conservation of sand grains, transfer 
of fixed number of grains at each toppling and isotropy. 
Thus, the results would be equally applicable to models 
in which toppling conditions are different, or the transfer 
of particles is stochastic, as in the Manna model. 

In Fig 3., we have shown the results of a MC simu- 
lation study of the DRT in four different models: (a) 
the 1-dimensional BTW model for L=100, (b) the 2- 
dimensional BTW model defined on a cylinder of size 
50x50, (c) the 1-dimensional Manna model with i = 100 
with rule that if Zi exceeds 1, then 2 particles are trans- 
ferred, each randomly to one of the neighboring sites, and 
(d) the 2-dimensional Manna model on 50 x 50 cylinder 
with two grains transferred at each toppling, each grain 
in randomly chosen direction. The number of particles 
used in each simulation was 10*^. We plot the probabil- 
ity of survival of the marked grain as a function of time 
t/M, where M was determined from the simulation di- 
rectly. We see a very good collapse and agreement with 
the theoretical prediction that exp{—t/M) for i 3> 1. 

For small t, the probability distribution is determined 
by the grains which are added very near the boundary. 
Boundary avalanches are not properly taken care of by 
our analysis. In particular, it is easy to see that the prob- 
ability of added grain coming out immediatly is nonzero 
in d-dimensions, and varies as l/L for large L (this be- 
ing the ratio of surface to volume). But Ea. (|13|) would 
give this to be 0{L~'^). This comes from the fact that 
near the boundary the height distribution is modified, 
resulting in the effective K becoming different near the 
boundary. Ea. H13(l is valid for i ^ 1. 

In Fig 4., we have plotted the difference 5S{t\L) be- 
tween the survival probability from MC simulation data 
and the scaling theory prediction [Eq. (|14|l ]. for a 2- 
dimcnsional BTW and the Manna models on an L x L 
square lattice for different values of L. We find that 
the curves for different L collapse onto each other if 
L5S{t\L) is plotted versus t/Lp' , indicating that the cor- 
rection 5S{t\L) to the Eq. H14() has the scaling form 
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5S{t\L) ^ -git/L') 



(15) 



P{x,t) = ^exp{-Kt/2V). 



(13) 



where the correction-to-scaling function .9(2;) is different 
for the two models. 

For the one-dimensional BTW model, with grains 
added everywhere with equal probability, we can de- 
termine exactly the leading 0{1/L) correction to the 
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scaling solutions Ea. (|14() . Thus, from the scaling solu- 
tion Eq.lO, we get Prob{T = 0\L) = 1/L + 0{L-^). 
But a straightforward calculation shows that actually 
Prob{T = 0\L) = 2/L + 0{L-^). However, Prob{t\L) 
for t = 1,2,..., is correctly given to the lowest order 
by 1/L. Assuming that the remaining distribution is a 
simple exponential, we get 

Sid BTw{t) = ^<5*^o + (1 - \)exp[-^{l - |)], (16) 

where we have used the normalization condition S'(O) = 
1, and added slO{1/L) correction term to the coefficient 
in the exponential. Using the condition that first moment 
of this distribution is M = + 1), we get a — 1/2. 

In Fig. 5, have shown the results of a Monte Carlo simu- 
lation of this case for L = 40 and number of grains= 10^, 
and compared with the theoretical prediction [Eg. 1)16(1 ]. 
We see that the agreement is very good. 

It is interesting to note that the DRT [Ea. H13() ] has a 
simple universal form, and does not depend on the critical 
exponents for the distribution of avalanches, which differ 
for BTW and Manna models, and also depend on dimen- 
sion in a nontrivial way. In contrast for the Oslo ricepile 
model, Christensen et al and Boguna and Coral |l7| 
found that the DRT for the 1-dimension ricepile of size 
L, with L 3> 1, does involve nontrivial exponents, and 
has the form 

Proh{T\L) = j^f{T/U). (17) 

The exponent ~ 1.3 is related to the roughness of the 
ricepile surface. The function f{x) takes a constant value 
for X small, and varies as x~^ for large x, with b ~ 2.4. 

To study such models theoretically, we have to consider 
sandpile models where there is a stack at each site into 
which grains are put in, and there are specific rules about 
which grains leave the stack. For example, one could as- 
sume that incoming grains are put on the top of the stack, 
and outgoing ones leave from the top. This corresponds 
to the last-in-first-out rule. Alternatively, we can choose 
the first-in-first-out rule, or choose the outgoing grains 
probabilistically, with probability of selection depending 
on distance of grain from top of the stack. Clearly, the 
DRT will have different behavior for different rules. We 
hope to address this question in the future publication. 

For the Oslo ricepile model, if we consider only grains 
those are not permanently stuck in the pile as consti- 
tuting "active mass" of the ricepile, all mass above the 
minimum slope of the pile is active. The configuration 
with the minimum slope is recurrent, and will recur in- 
finitely often in the steady state as grains are added to 



the pile. Therefore, all the grains added after the pile has 
reached the minimum slope have a finite residence-time 
in the pile. Then the argument given earlier in this paper 
implies that mean active mass of the ricepile should vary 
as i^. But the result, obtained in and 0|, shows 
that it varies as L^-^ . The reason for the discrepancy in 
the estimate of mean mass in 0, 0| is presumably due 
to the very long residence of the grains which happen to 
get deeply embedded, making the estimate of the first 
moment of the DRT unreliable from short simulations. 



We thank A. Nagar for discussions. 
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FIG. 1: Motion of three grains starting from x = 20, 50, 80 

on a one dimensional sandpilo of length L = 100 where sand 
grains are added only at the right end. 
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FIG. 2: Survival probability versus residence time t for the 1- 
dimensional BTW model in two cases with grains added only 
at one or both ends. The theoretical result (full curve) and 
the simulation result (dotted line) match perfectly. L = 100 
for both the cases. 
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FIG. 3: Semilog plot of the survival probability of the marked 

grain as a function of scaled residence time t/M for four dif- 
ferent cases; (1) one dimensoinal BTW model, (2). one di- 
mensoinal Manna model, (3) two dimensoinal BTW model, 
(4) two dimensoinal Manna model. We have chosen L = 100 
in one dimension and 70 x 70 cylindrical square lattice in two 
dimensions. 
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FIG. 4: The scaled correction term L5S{t\L) versus t/L^ in 
2-dimensional BTW and Manna model from simulation for 
three different values of L = 13, 20, 30. The curve with higher 
peak is for the Manna model. 
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FIG. 5: Finite size correction for the one-dimensional BTW 
model with grains added everywhere for L — 40. Plotted 
is the scaled deviation, L5S{t/L), from the simple exponen- 
tial scaling solution [Eg. 1141 1. of the simulation data and the 
theory [Ea. iT?H ]. 



